Space charge limited flow of a thin electron beam confined by a strong magnetic field 
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An approximate analytic theory is developed and implemented numerically for calculating the space 
charge limited current and electric field of a thin cylindrical beam or current sheet between two wide 
parallel electrodes. The flow is confined by a sufficiently strong magnetic field. Assuming that the 
potential and current density are almost homogeneous in the direction transversal to the flow we 
compute the beam current and profile by a variational method. We find that the average current 
density scales as the reciprocal of the beam width when the latter becomes very small. The total 
cylindrical beam current thus decreases proportionly to its diameter while the total current of a 
sheet becomes almost independent of the width in this regime. 

PACS: 52.27.Jt; 52.59.Sa; 52.59.Wd; 85.45.Bz 
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1. Introduction 

Space charge limited electron flow in two and three 
dimensions presents a difficult nonlinear mathematical 
problem whose solution is important for many practical 
applications. In the design of high power electron beams 
[1-4] the one dimensional Child-Langmuir limit (CLL) 
[5] has been a benchmark for almost a century but cor- 
rections have to be made for the "current wings" near 
the boundaries of the flow. Consequently the current of 
narrow beams [6-8] show great divergences from CLL. 

In a previous article [9] we considered a planar emitting 
region whose width 2a is much larger than the cathode- 
anode distance (1 in our dimensionless units). In the case 
of narrow beams, 2a <C 1, considered here, the moderate 
variations of the potential in the transversal direction will 
be used to reduce the dimensionality of the problem. We 
assume that the two parallel electrodes arc large enough 
for the electric field to be homogeneous far away from 
the beam. A magnetic field directed along the y-axis 
keeps the current, which is emitted either by a long nar- 
row straight strip or a circular disk, perpendicular to the 
electrodes without spreading out as shown in Fig. 1. 

The outline of the paper is as follows. In section 2 we 
set up the problem and solve the Laplace equation for 
the vacuum field outside the current sheet. This solution 
serves as a boundary condition to derive a closed equation 
for the potential and current in the space charge region 
using the continuity of the electric field. In section 3 we 
specialize to the case of narrow beams and make certain 
approximations. We then solve the approximate problem 
iteratively in section 4 using a direct variational method. 
The results are presented in section 5. We also consider 
there the asymptotics of the current density in the limit 
a — > 0. In section 6 we extend our approach to cylindrical 
current beams. Section 7 is devoted to a discussion of 
results. 



2. Formulation of problem. 

We will use dimensionless units: linear sizes are mea- 
sured as fractions of the inter-electrode distance D, the 
potential (j>(x, y) is given as a fraction of the inter- 
electrode voltage V, and the current density j(x) is in 
the units of the Child-Langmuir current density jo — 
V 3 / 2 V2e/9TTD 2 y/m, obtained in the limit a = oo [5]; e 
and m are the electron charge and mass. 
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FIG. 1. Cross-section of the current with an expanded 
schematic graph of a contour of constant potential and electric 
field 

The nonlinear Poisson equation which determines both 
the potential <f>(x, y) between the electrodes and the cur- 
rent density j(x) for a current sheet, see Fig. 1, has the 
form [9] 
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d 2 



Aj{x) 



dx 2 dy 2 9^0^) 



(1) 



Eq. (1) is to be solved with the boundary conditions 

<p(x,0) = 0, 4>(x,l) = I. (2a) 

The unknown space charge limited current j(x), in the 
region \x\ < a, is to be determined from the condition 
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dy 



(x,0) = 0, for \x\ < a. 



(26) 
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For |x| > a, j(x) — and (1) reduces to the Laplace 
equation in vacuum. We assume that (1) and (2) are 
well posed but leave a mathematical investigation of this 
point for the future. 

A natural approach to solving (1) is to separate it 
into a standard boundary value problem for the Laplace 
equation outside the space charge region, \x\ > a, and 
a nonlinear inner problem in the rectangle, —a < x < 
a , < y < 1. Taking the potential ipiy) = <?K a jZ/) 
on the sheet surface as given and using the symmetry 
4>{x, y) = (j){—x, y) wc find <f>(x, y) for x > a. From this 
one calculates the external electric field in the x-direction 
e+{y) = — ^(a,y) at x = a + . The continuity of the po- 
tential and electric field at the space charge boundary 
gives the additional boundary conditions tp(y) = </>(a,y) 
and e_(y) = e + (y) — e(y) needed for solving (1) in the 
rectangle \x\ < a, < y < 1. 

Vacuum Fields 

We want to solve the Laplace equation (1) in the region 
x > a with the boundary conditions, 

<f>(x, 0) = 0, 4>(x, 1) = 1, 0(a, y) = ^(y), (3) 



and 



(x, y) — > y for x — > oo, 



expecting ip(y) to be monotone with finite first deriva- 
tives. Taking the potential in the form of a Fourier series 

_ oc 

<f>(x, y) = y + - V -V fe («-z) sin nky, x>a (4) 



k=l 



(3) is satisfied if 
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(y) cos irkydy. 



(5) 



An illustrative (and useful for our computations) ex- 
ample is the case ip(y) = J/ 4 ^ 3 - Eq. (5) then gives 



(—ikn) 4//3 7 ( -, — ikir 



where j(u, v) is the incomplete gamma function [10]. 

Using (4) and (5) the x-component of the electric field 
at the beam boundary x — a can be written as 



e(y) = _ ^( a '2/) = 2^a fe sin7rfcy. 



(6) 



k=i 



Substituting (5) into (6) and summing up the series ex- 
plicitly we get an alternative form for e(y) 



e(y) = sin7ry 



I 



1 i>'(t)-ny) 

cos nt — cos ny 



dt. 



(7) 



We note that Eq.(7) can also be derived by using confor- 
mal mapping, as in [9]. For ijj(y) = y 2 Eq.(7) becomes 



4 f y 

<y) = - \ 

7T Jo 



mtan -^dt. 



The inner reyion 

Eq. (1) together with boundary conditions (2) and (7) 
now define the space charge limited current problem in 
the rectangle \x\ < a, < y < 1, the shaded region in Fig. 
1. Also shown there is an expanded cross section together 
with a schematic drawing of an equipotcntial surface and 
electric field vector with its x-componcnt. As already 
noted the CLL potential (a = oo) is 4>q(x, y) = y 4 ^ 3 and 
j(x) = 1. For finite width of the current sheet <f>(x, y) — ► 
<pi(y) = y as \x\ — ► oo, and we expect <j){x,y) inside the 
beam, \x\ < a, to satisfy the inequality y 4 / 3 < (f>(x,y) < 

y- 

3. Narrow Beams 

Integrating (1) between x — and x — a, we obtain 



Q2 ra 4 ra 

Eq.(8) can be rewritten as 



-e(y) + a- 



(xi,y) ^aj{x2) 



dy 2 



dx. (8) 



(9) 



where < x\(y),X2{y) < a. Assuming now that for 
a <C 1 there is only a small variation in 4>(x, y) for 
< x < a, at least when y is not too close to zero or one, 
we approximate (9) by setting <p(xi,y) = <p(x2,y) — 4>{y) 
and j(x2) = j, where j is a parameter, close to the av- 
erage current density, which will be determined together 
with (f>(x,y). We expect that in the limit a — > 0, e(y), 
ip{y) and aj will become independent of a (see section 
6). We therefore take the equation 



-e(y) + a 



d 2 ^y) 

dy 2 



4aj 



(10) 



together with (7) as a suitable approximation for »« 1 
and try to solve them numerically. The boundary condi- 
tions for tp come from (2): ip(0) = 0, ip'(0) = 0, ip(l) = 1. 

An alternate way to arrive at (10), which also gives 
some information about the beam profile, is to approxi- 
mate (f)(x, y) for a <C 1 and y not too close to zero or one 
, as 

<f>{x,y) c ^y) + ^l-^ e (y) (11) 
so that e(x, y) ~ | e(y). In the same spirit we would have 

^2 s 



j{x) ~ j[r] + i 



(12) 
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where j(f? + £/3) is the mean current density and £ repre- 
sents the relative height of the current wings (expecting 
X] to be close to 1). Substituting (11), and (12) in (1) we 
obtain two ordinary equations, namely (10) and 
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i 2 — x 2 d 2 e 



rj + £x 2 /a 2 



8aj dy 2 ^ + e{y)(a 2 - x 2 )/2a V? 
The ratio 



-TT- ( 13 ) 



C(y) 



(14) 



gives the relative depth of the potential profile (11) at its 
deepest point (x = 0) as a function of y. 

4. Numerical Algorithm 

To solve (10) numerically we write ip in the form 



*P(y) = y i/3 + E c ™/«^)' 



(15) 



where the f m are a set of twice differentiable functions 
with the properties / m (0) = /„(0) = f m (l) = 0, m = 
1, 2, n to satisfy the boundary condition (2) for ip(y). 
The boundary electric field is similarly expressed as 

n 

e(y) = e (y) + ^ (y)- (16) 

m— 1 

For computing e m (y), m = 0, 1, .., n, we use (6) with the 
coefficients, given by (5), 



Jo 



cos nkydy. 



(17) 



After this we apply the direct variational method to min- 
imize the functional 3>(ci, C2, c„; j) 



e(y) - a 



d 2 ip(y) 4aj 



y 4/3 A/ (18) 



in terms of the parameters c m and j. The factor y 4 / 3 is 
used in (18) for regularization. "Neutralization" of the 
nonlinearity of l/y / tp(y) in (18) is achieved by an iter- 
ation procedure where at each step we use parameters 
ci,C2,...,c„ evaluated in the previous step. The func- 
tional $ thus becomes bilinear in terms of its parameters 
and we solve repeatedly a linear algebraic system 



= o, 



dc m 



= 0, m = 1, 2, .., n 



(19) 



until its solutions stabilize. 

A serious problem in this variational computation is a 
good choice of the trial functions f m . We want in partic- 
ular the right behavior near the singular points y = and 
y = 1. It is easy to see from (6) that e(0) = e(l) = 0, but 



their derivatives are generally divergent. On the other 
hand if the error in a small neighborhood of these points 
is not large their overall impact on e(y), < y < 1 is not 
significant. To eliminate these regions in computing $ we 
set y m in ~ 0.09 and y max w 0.99 in (18). We monitor the 
validity of our algorithm by two indicators: (1) the mini- 
mum value of $ after the iterations get stabilized and (2) 
by deviations of the left part of Eq. (10) multiplied by 
9-\/0(y)/4afrom the current j, at y = 0.1, 0.2, 0.9. For 
flexibility in approximating tl>(y) and in order to keep the 
matrix of the system (19) away from any degeneracy we 
choose for f m (y), m = 2,3, 10 a set of partially over- 
lapping functions of finite support. The function f\(y) 
which corrects the CLL fo(y) for small y is taken of the 
form e-Pvy 4 / 3 with (3 ~ 40 - 60. 

To carry out our scheme practically we calculate in ad- 
vance the partial electric fields e m (y) on a homogeneous 
grid of 1000 points and then apply the iterations. These 
usually converge rapidly (after less than 10 steps) while 
each iteration takes a very short time. Some series con- 
verge rather slowly (as k~ 4 ^ 3 ), but even for a precision of 
10~ 9 , i.e about 10 5 terms in the series, the computation 
of all partial fields is very fast. 

After finding ip(y) and e(y) the same, but much sim- 
pler, procedure can be applied for evaluating the parame- 
ters rj and £ in (12). We use for Eq.(13) the least squares 
method, where the weight function in the functional is 
chosen in such a way as to eliminate the derivatives of 
e(y), and the computation is done around the middle of 
the beam. 

5. Results for the current sheet. 

In Table 1 we present the results of the variational 
procedure described in the last section for different values 
of a. We show there the parameter j which is found 
by (19) and the mean deviation A of the corresponding 
quantity, defined in terms of ip(y) and e(y) in Eq.(10), on 
the segment 0.1 < y < 0.9. The precision of computation 
is controlled by the ratio $/$ at the stationary point, 
where 



f)7 



4/3 



i>(y) 



dy. 



a 


0.2 


0.1 


0.05 


0.02 


0.01 





j *a 


0.549 


0.447 


0.394 


0.360 


0.347 


0.326 


A 


0.074 


0.080 


0.075 


0.062 


0.053 


0.041 




.0052 


.0061 


.0055 


.0043 


.0044 


.0067 




0.042 


0.032 


0.021 


0.009 


0.004 






0.026 


0.021 


0.014 


0.006 


0.003 






TABLE 1. 


Flat current sheet 





Also shown there are the depth £ of the current density 
dip across the beam and the quantity 



"-3' 



3 



which gives the relative difference between the value j 
computed by solving Eq.(lO) and the mean current den- 
sity j(x). The last column in Table 1 exhibits the limiting 
value of these quantities when a — > using a procedure 
described below. 

The limit a — > 0. 

Our model assumes that the Larmor radii of electrons 
do not exceed a and this makes the limit a — > physi- 
cally questionable since this would require the magnetic 
field to become unreasonably strong. Nevertheless the 
limit a — > is very interesting both mathematically and 
practically for evaluation of the prefactor in the scaling 
law for j. 

Looking at Eq. (9) we expect that in the limit a — > 
the term aj(x2) — * A, independent of a, ip(x2, y) — > il>(y), 



and 



. d 4>(x,y) 



except at y = 0, y = 1. The limiting 



ip(y) and e(y) will then satisfy the equation 



e(y) + 



4A 



(20) 



with boundary conditions which we do not know. We 
therefore solved (20) in a truncated interval y\ — 0.09 < 
y < 0.998 = 2/2 using the same routine as before. The 
result of the computation is presented in the last column 
of Table 1. We believe now that by improving and ex- 
tending the set of trial functions f n {y) one can gradually 
take j/i, 2/2 closer to and 1 respectively, but we think 
that A ~ 0.3 is a good approximation. 

To get an "exact" equation for tp and A we integrate 
Eq. (20) over an interval [yi , 2/2] < 2/1 < 2/2 < 1 and sum 
the series in (5) and (6) to yield a nonlinear eigenvalue 
problem 
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47T 



ip'(t) In 



cos 717/1 — cosnt 



COS 7T2/2 — cos irt 



dt 



f y 

Jyi 



dt 



(21) 

where the boundary conditions are not needed. We post- 
pone investigation of (21) for the future. 

6. Cylindrical beam 

We consider now the more common case of electron 
beams of a compact cross section when the emitting part 
of the cathode as well as the cross section of the beam 
are circular. In cylindrical coordinates Eq. 1 takes the 
form 



ld_ 

r dr 



dr 



d 2 <j> 

dy 2 



4j(r) 



^<j>(r,yY 



(22) 



with the boundary conditions 



0(r, 0) = 0, 0(r, 1) = 1; ^(r, 0) = 0, for r < a (23) 

and j(r) = when r > a. Fig.l represents now the beam 
cross section which passes through its axis r = 0. 



Carrying out again a Fourier expansion for <p(r, y) in 
the vacuum region r > a yields the potential as a series 

^, y )=2/+-g T ^ R sm7rfc2/,r> fl (24) 

where can be found by (5) with ip(y) = (j)(r,y) at 
r = a. Calling now e(y) the radial component of the 
electric field, at r = a, leads to the analogue of (6) 



«(») = -f(«.») = 2E 



Ki(ytka) 
^ K (Trka) 



aksinnky. (25) 



Here and in (24) Kq, K\ are the modified Besscl functions 
which decay exponentially at infinity. 

We can apply the same technique as in (9)- (12) and 
obtain instead of Eq.(10) 



-2e(2/) +a ^ 



Aja 



(26) 



Eqs.(13) and (19) for the current and potential profiles 
are the same. The factor 2 in (26) shows the stronger 
effect of the surrounding electric fields on the cylindrical 
beam compared with that on the narrow sheet. The same 
numerical scheme as in part 4 yields the results shown in 
Table 2 



J * a 
A 



0.2 
5.574 
0.090 
.0051 
.1089 



0.1 
5.532 
0.069 
.0029 
.0566 



0.05 
5.452 
0.052 
.0023 
.0293 



0.02 

5.422 

0.037 

.0026 

.0119 



0.01 

5.410 

0.036 

.0029 

.0059 



TABLE 2. Cylindrical beam 



The current density behaves like j <~ 5.4/a and there- 
fore the total beam current becomes proportional to a, 
when a« 1. The prefactor for cylindrical beams is sig- 
nificantly higher than the one for current sheets because 
in this case the electron flow is completely surrounded 
by the vacuum field, but a compact expression like (21) 
for the limiting parameter A and the boundary potential 
ijj is impossible here. The profiles of the current density 
and potential are flatter than in the current sheet. 

7. Discussion 

As already noted the current density j scales approx- 
imately as a -1 and thus will be large when the beam 
is narrow with the total current of the sheet becoming 
independent of the sheet thickness (provided by the un- 
limited emissivity of the cathode). The current density 
grows slightly near the beam boundaries, but even for 
the flat sheet whose thickness is 0.4, which is not small, 
this rise is less than 10% (the evaluation of this quantity 
is more reliable for a < 0.1). On the other hand when 



4 



y is small the potential across the beam width varies 
much stronger. In Fig. 2 we show the electric field e(y) 
for a = 0.1 and the maximum deviation £(y) (14) of the 
cquipotential surface from the horizontal for a — 0.1 and 
a = 0.02. 

When a — > the potential ip(y) approaches the so- 
lution of (21). This function does not differ too much 
from y 4 / 3 . Only when y is small does ip(y) increase sub- 
stantially. It has also an irregular behavior near y = 1 
where ip"(y) ~ 4j/9, since e(l) = 0. On the other 
hand ip"(y) is of order of 1 almost everywhere on the 
segment (0, 1) because if tp"(y) > C, C > on an inter- 
val y\ < y < i/2 then, using the inequality ip'(y) > 0, we 
would get VK2/2) — V'(yi) > C(2/2 — 2/i) 2 /2 which, if C is 
large, will contradict the condition ip(y) < 1. Therefore 
the term e(y) is dominant almost everywhere in the left 
side of (10) for a < 1. 
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FIG. 2. Current sheet. Plots of C(y) for a = 0.1, 0.02 (solid 
lines) and e(y) for a = 0.1 (dashed line) 

The cquipotential surfaces deviate from horizontal 
planes only when y is small (in the case of a = 0.01 
this effect is of lower order). Clearly our method of com- 
putation, which neglects the x-derivatives higher than 2, 
is not accurate near y — where the electron density has 
a complicated structure in the x-direction. We believe 
nevertheless (see discussion in part 6) that these short 
intervals do not play a crucial role because the current 
is limited by the whole space charge distributed between 
the electrodes. 

There is a simple generalization of the parabolic shape 
assumption used in (11), (12) for <f> and j{x) for larger o 
[11]. This uses for the approximate solution of Eq.(l) the 
following ansatz 

<j>(x, y) = ip(y) + e(y)w(x), j(x) = j[r) - £w(x)], (27) 

where w(x) = ( cosh ga — cosh gx)/ksinh ga and the posi- 
tive parameter g can be taken from our previous work [9] 
(g w 3.88). Substituting (27) into (1) yields two ordinary 



differential equations similar to (10) and (13). The ad- 
vantage of the scheme used here for narrow beams is the 
transparent relationship between the exact Eq.(9) and 
our main Eq.(10) which is approximate for a ^ but 
becomes exact, see (20), when a — > 0. 
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